##########################################################################################
#### 20231211 ArchR分析步骤 ####

options(stringsAsFactors=F)
library(Seurat)
library(ArchR)
library(optparse)
library(BSgenome.Hsapiens.UCSC.hg38)

##########################################################################################

option_list <- list(
    make_option(c("--rna_data_file"), type = "character"),
    make_option(c("--input_dir"), type = "character"),
    make_option(c("--cluster_cell_file"), type = "character"),
    make_option(c("--cpu"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
  
  cluster <- "cluster0"
  cell_type_file <- "~/20231121_singleMuti/config/cluster_celltype.csv"
  input_file <- "~/20231121_singleMuti/results/atac_res/testis_combined_peak.Rdata"
  out_path <- "~/20231121_singleMuti/results/atac_res"

}

dir.create(out_path , recursive = T)
dir.create(paste0( out_path , "/motif_images/testis_combined_peak.motif.Rdata" ) , recursive = T)
setwd(out_path)

##########################################################################################

a <- load(input_file)
cell_type <- fread(cell_type_file , header = F)

##########################################################################################

use_cell <- subset( cell_type , V1==cluster )$V2

projHeme4 <- testis_combined_peak
projHeme4@peakSet <- projHeme4@peakSet[which(names(projHeme4@peakSet) == use_cell),]

projHeme5 <- testis_combined_peak

##########################################################################################
## https://www.archrproject.com/bookdown/motif-deviations.html
projHeme5 <- addMotifAnnotations(ArchRProj = projHeme5, motifSet = "cisbp", name = "Motif")

## 我们还需要添加一组用于计算偏差的背景峰值。使用 chromVAR::getBackgroundPeaks() 函数选择背景峰，
## 该函数使用马哈拉诺比斯距离根据所有样本中 GC 含量和片段数量的相似性对峰进行采样。
projHeme5 <- addBgdPeaks(projHeme5)

## 计算所有主题注释中的每个细胞偏差
projHeme5 <- addDeviationsMatrix(
  ArchRProj = projHeme5, 
  peakAnnotation = "Motif",
  force = TRUE
)

out_file <- paste0( out_path , "/testis_combined_peak.addmotif.Rdata" )
save( projHeme5 , file = out_file)

## 等矩阵中，我们将染色体信息存储在 . 没有任何相应的位置信息，
## 而是使用两个不同的序列名称将 chromVAR 中的“偏差”和“z 分数”存储到同一矩阵中
#plotVarDev <- getVarDeviations(projHeme5, name = "MotifMatrix", plot = TRUE)

# Get locations of motifs of interest:
motifPositions <- getPositions(projHeme5, name="Motif")
motifGR <- stack(motifPositions, index.var="motifName")

##########################################################################################

if(1!=1){
  motifs <- c("GATA1", "CEBPA", "EBF1", "IRF4", "TBX21", "PAX5")
  markerMotifs <- getFeatures(projHeme5, select = paste(motifs, collapse="|"), useMatrix = "MotifMatrix")

  p <- plotEmbedding(
      ArchRProj = projHeme5, 
      colorBy = "MotifMatrix", 
      name = sort(markerMotifs), 
      embedding = "UMAP",
      imputeWeights = getImputeWeights(projHeme5)
  )

  p2 <- lapply(p, function(x){
      x + guides(color = FALSE, fill = FALSE) + 
      theme_ArchR(baseSize = 6.5) +
      theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
      theme(
          axis.text.x=element_blank(), 
          axis.ticks.x=element_blank(), 
          axis.text.y=element_blank(), 
          axis.ticks.y=element_blank()
      )
  })
  p3 <- do.call(cowplot::plot_grid, c(list(ncol = 3),p2))

  out_file <- paste0( out_path , "/motif_images/testis_combined_peak.motif.pdf" )
  ggsave(out_file , p3)
}